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ABSTRACT 

We present a multi-dimensional numerical code to solve isothermal 
magnetohydrodynamic (IMHD) equations for use in modeling astrophysical 
flows. First, we have built a one-dimensional code which is based on an explicit 
finite-difference method on an Eulerian grid, called the total variation diminishing 
(TVD) scheme. The TVD scheme is a second-order-accurate extension of the Roe-type 
upwind scheme. Recipes for building the one-dimensional IMHD code, including the 
normalized right and left eigenvectors of the IMHD Jacobian matrix, are presented. 
Then, we have extended the one-dimensional code to a multi-dimensional IMHD code 
through a Strang-type dimensional splitting. In the multi-dimensional code, an explicit 
cleaning step has been included to eliminate non-zero V • B at every time step. 

To test the code, IMHD shock tube problems, which encompass all the physical 
IMHD structures, have been constructed. One-dimensional and two-dimensional 
shock tube tests have shown that the code captures all the structures correctly 
without producing noticeable oscillations. Strong shocks are resolved sharply, but 
weaker shocks spread more. Numerical dissipation (viscosity and resistivity) has been 
estimated through the decay test of a two-dimensional Alfven wave. It has been found 
to be slightly smaller than that of the adiabatic magnetohydrodynamic code based on 
the same scheme. As an example of astrophysical applications, we have simulated the 
nonlinear evolution of the two-dimensional Parker instability under a uniform gravity. 
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1. INTRODUCTION 

Over the last two decades, conservative upwind differencing schemes have proven to be very 
efficient for solving adiabatic hydrodynamic and magnetohydrodynamic (MHD) equations. These 
methods generally depend on the calculated estimates of mass, momentum and energy fluxes as 
well as magnetic field flux across cell boundaries based on the so-called "Riemann" solutions from 
the basic conservation laws. Examples for hydrodynamics include extensions of Godunov's scheme 
(Godunov 1959), such as the MUSCL scheme (Van Leer 1979) and the PPM scheme (Colella & 
Woodward 1984), as well as those based on approximate flow eigenstates such as the Roe's scheme 
(Roe 1981), the TVD scheme (Harten 1983) and the ENO scheme (Harten et al. 1987). Works for 
MHD include Brio & Wu (1988), Zachary & Colella (1992), Zachary, Malagoli, & Colella (1994), 
Dai & Woodward (1994a,1994b), Ryu & Jones (1995) (RJ, hereafter), Ryu, Jones, & Frank (1995) 
(RJF, hereafter), Powell et al. (1995), and Roe & Balsara (1996). Brio & Wu applied the Roe's 
approach to the MHD equations. Zachary and collaborators used the BCT scheme to estimate 
fluxes for the MHD conservation equations. Dai &: Woodward adapted the PPM scheme to MHD. 
Ryu and collaborators extended the Harten's TVD scheme to MHD. Powell and collaborators 
developed a Roe-type Riemann solver with an eight-wave structure for MHD. Roe & Balsara 
constructed one variety of linearized Riemann solutions for MHD. The upwind schemes generally 
share an ability to sharply and cleanly define fluid discontinuities, especially shocks, and exhibit a 
robustness that makes them broadly applicable. 

The assumption of adiabatic flows holds in the limit where cooling is negligible or the cooling 
time scale is much larger than the dynamical time scale. However, in the other limit where the 
cooling time scale is much shorter than the dynamical time scale, the assumption of isothermal 
flows becomes physically more valid (e.g., Draine & McKee 1993 and references therein). Of 
course, if cooling time scale is comparable to dynamical time scale, cooling should be considered 
explicitly. 

Usually, numerical simulations of isothermal flows are made with adiabatic codes by setting 
the adiabatic index, 7, close to unity. Truelove et al. (1998) showed that with 7 as close as to 
unity as 1.001, their adiabatic hydrodynamic code can follow isothermal collapse without any 
significant deterioration of accuracy. We also observed that with 7 = 1.001 the adiabatic TVD 
MHD code (RJ) captures structures in isothermal magnetohydrodynamic (IMHD) shock tubes 
without noticeable error. Yet, it is desirable to build codes specifically for isothermal flows, since 
those codes are simpler and faster than adiabatic ones. That is because the energy conservation 
equation need not to be solved in isothermal codes. As a result, the entropy mode, which carries 
the contact discontinuity, need not to be considered. In the current paper we describe an IMHD 
code based on Harten's TVD scheme. It is the same scheme that was used for the adiabatic MHD 
code in RJ and RJF. Balsara (1998b) developed an IMHD code also based on an upwind scheme, 
but his scheme is different from ours. 

In §2, we give recipes for the development of one- and multi-dimensional IMHD codes. In 
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§3, we present the results of tests that include one-dimensional and two-dimensional shock tube 
problems, the decay of an Alfven wave, and the nonlinear evolution of the Parker instability under 
a uniform gravity. Conclusions follow in §4. 



2. NUMERICAL SCHEME 



2.1. The Equations for Isothermal Magnetohydrodynamics 



MHDs describes the behavior of the combined system of a conducting fluid and magnetic 
fields in the limit that the displacement current and the separation between ions and electrons 
are neglected. So the MHD equations represent coupling of the equations of fluid dynamics with 
Maxwell's equations of electrodynamics. By ignoring the effects of electrical resistivity, viscosity, 
and thermal conductivity, and imposing isothermality on the conducting fluid, we get the following 
IMHD equations: 

dp 



— + v ■ Vv + -V(a 2 p) 
at 



P 



1 
P 



0, 



(VxB)x8 = 0, 



dB 
~dt 



with an additional constraint 



V x (v x B) 



V -B = 0, 



(2-1) 
(2-2) 
(2-3) 

(2-4) 



for the absence of magnetic monopoles. Here a is an isothermal sound speed, and other notations 
have their usual meanings. We incorporate a factor of 1/yAjv into the definition of B so that the 
factor of 4-7T does not appear in Eq. (2-2). 



In Cartesian coordinates, the above equations are written in a conservative form as 

dq , 8F X , dF v , 8F Z n 



Ox 
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(2-5) 
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(2-6) 



/ 



with F y and F z obtained by properly permuting indices. With the state vector q and the flux 
functions F x (q), F y (q), and F z (q), the Jacobian matrices, A x (q) = dF x /dq, A y (q) = dF y /dq, and 
A z {q) = dF z /dq are formed. A system is called hyperbolic if all the eigenvalues of the Jacobian 
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matrices are real and distinct and the corresponding set of right eigenvectors is complete (Jeffrey 
& Taniuti 1964). The system of the ideal, adiabatic MHD equations is known as non- strictly 
hyperbolic, since some eigenvalues coincide at some points (Brio & Wu 1988; Roe & Balsara 1996). 
The eigen-structure of the IMHD equations, which is presented in the next subsection, is very 
similar to that of the adiabatic ones. It is easy to show that the IMHD equations also form a 
non-strictly hyperbolic system. 



2.2. One-Dimensional Code 



Our strategy for developing a one-dimensional IMHD code is based the TVD scheme (Harten 
1983) which was devised to improve the first-order-accurate Roe's upwind scheme (Roe 1981) 
into a second-order-accurate one. For it, we derive the eigenvalues and eigenvectors of the system 
of the IMHD equations, which are given below. With the eigenvalues and eigenvectors, it is 
straightforward to apply the construction procedure for the one-dimensional adiabatic MHD code 
(RJ) to an isothermal analogue. Even though the procedure is described in RJ, we here repeat it 
to make this paper self-contained. Special attention is given to the orthonormal eigenvectors of 
the system of the IMHD equations. 

We consider, as an example, plane-symmetric, one-dimensional flows exhibiting variation along 
the x-direction. Then y- and z-derivatives in Eq. ( |2-5[ ) are zero, and we have the one-dimensional 
IMHD equations 

dq 8F X 



+ 



0. 



— B 

dx x 



is %B a 



dt dx 

The fifth equation in the system of Eqs. (|2^ 

0. These imply that initially B x should be spatially 
constant and be kept constant during the evolution of the flow. So we need not include the 
equation for B x in a one-dimensional code. 



where q and F x are defined in Eq. 
and the constraint in Eq. (|2-4|) is 



(2-7) 
= 0, 



The Jacobian matrix dF x /dq of the system of Eqs. fl2-7| ) is given by 
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where b x ^ y ^ z = B x ^^ z j ' J~p. The six eigenvalues in non-decreasing order are 



(2-8) 



a i 



(2-9) 



0,2 = V x - C a , 



(2-10) 
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03 = V x - C s , 


(2-H) 


a 4 = v x + c s , 


(2-12) 


a,5 = v x + Co, 


(2-13) 


a 6 = v x + c f , 


(2-14) 



where c/, c a , c s are the fast, Alfven, and slow characteristic speeds, respectively. There is no 
entropy mode for the IMHD equations. The quantities ai, - ■ ■ ,oq represent the six speeds with 
which information is propagated locally by three MHD wave families. The three characteristic 
speeds are expressed as 



Cf = {\[a 2 + bl + bl + b 2 z + + hi + b\ + ft*)* - Aa%l\ } 

Ca = \b x \i 



C, = 



1 r 

2 



a 2 + K + K 



a 2 + bl + 62 + & 2)2 _ 4a 2 6 2 



} 



1/2 



1/2 



The six right eigenvectors corresponding to the six eigenvalues are 

/ 1 \ 
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(2-15) 
(2-16) 
(2-17) 



(2-18) 



(2-19) 



(2-20) 
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Near the point where either b x = or b y = b z = 0, the above right eigenvectors are not well defined, 
with some elements becoming singular. By re-normalizing the eigenvectors, the singularities can 
be removed. The renormalized right eigenvectors are 



c f 



±Ca 



Of 

a f (v x ± c f ) 
afVy =F Ci s f3 y b x 
OifV z =F a s (3 z b x 

a s f3 y Cf 

CtsPzCf 
VP 

f \ 


T&sign(& x ) 
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\ _ hL 



\ 



) 



a s 
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CfVP 
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c fVP 



where a's and /3's are defined by 



a f 



hi 



b v 
b z 



fig + 6* 



At the points where b y = b z = 0, /3's are defined as a limiting value, i.e., 



(2-21) 



(2-22) 



(2-23) 



(2-24) 

(2-25) 
(2-26) 
(2-27) 

(2-28) 
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Similarly, at the point where b y = b z = and b 2 x = a 2 , a's are defined as 

af = a s = l. (2-29) 

The left eigenvectors, which are orthonormal to the right eigenvectors, Li ■ R m — Si m , are 

,(1) ,(2) ,(3) ,(4) ,(5) ,(6) 



r _ (iV L ) lK*) iW ]W /W \ 

1J V x ±Cf — \''V x ±Cfl''V x ±Cf1 l V x ±Cf1 l V x ±Cfl''V x ±Cf1 l V x ±Cf)l 

l i%c f = 7T a f a2 ± 7T [~OLfav x + a s c s (P y Vy + /? z v z )sign(6 :E ) 



n V2 



z( 2 ) I 1 



,(3) 1 



W± c/ = T7-a s /? ?/ c s sign(6 :E ), 



?2 
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T - (lW /( 2 ) /( 3 ) /( 4 ) 7( 5 ) 7(6) ^ 
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lv]±c a = ±\(PzV y - p y v z )sign(b x ), 



/ (2) _ () 
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C]±c s = ■T a s ( ?f T 7T [a s c a v x + a f Cf(P y Vy + /3 z w z )sign(6 :E )] , 

fl #2 

/ {2) - + — r> r 
*2±c = ±^a//?j/C/sign(6 ;E ), 



1 
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l { £ ±Ca = ±±-a f (3 z c f sign(b x ), (2-48) 

l ( 2±c s = ~a f fiyCfy/p, (2-49) 

l { Sc 3 = -ja f p z c fy /p, (2-50) 

where 

9 1 = 2(aja 2 + a 2 s c}), (2-51) 

9 2 = 2(a}c f a + a 2 s c a c s ). (2-52) 

Some elements in the normalized right and left eigenvectors are not continuous. In order to 
force them to be continuous, 

f 1, if b y >0 or b y = and b z > 
& V TJ \ -1, if 6 y < or ^ = and b z < 0, V ; 

is multiplied to R Vx ± Cs and L Vx ± Ca if a 2 > c 2 , and to R Vx ± Cf and L Vx ± Cf if a 2 <c 2 . 

Note that our eigenvectors have a different form from those derived in Balsara (1998a). The 
exact form of the eigenvectors does not matter, once all the singular points are taken care of. 

Here, we use the conventional indices. The superscript n represents the time step. The 
subscript i indicates quantities at the cell center, while i + \ marks those at the right-hand cell 
boundary. The subscript k represents the characteristic fields, with the order that k = 1 is for the 
field associated with eigenvalue v x — Cf, k = 2 for the field with v x — c ffi , k = 3 for the field with 
v x — c s , k = 4 for the field with v x + c s , k = 5 for the field with v x + c a , and finally k = 6 for the 
field with v x + Cf. 

An important step in the Roe's scheme (1981) is to determine a Roe matrix A x i+1 / 2 (qi, Qi+i) 
at the cell boundary from the adjacent state vectors, which satisfies the Roe's suggested properties. 
One of them is F Xj i + i — F x ^ = A X j +1 / 2 (qi+i — %)■ For the systems of the adiabatic and isothermal 
hydrodynamic equations, there exists a Roe matrix evaluated at the y/p-weighted average state 
(Roe 1981; LeVeque 1997). For the system of the adiabatic MHD equations, there is, however, 
no simple form of the Roe matrix except for the case with an adiabatic index 7 = 2 (Brio & 
Wu 1988). We have failed to find a simple form of the Roe matrix for the system of the IMHD 
equations, too. So we use an arithmetic averaging for the flow quantities at the cell boundary in 
the IMHD code, which was shown to work well in the adiabatic MHD code (RJ): 

Pi + Pi+i , n K .s 

Pi+1/2 = g ' ( 5 ) 

v x ,i+i/2 - 2 ' i^-ooj 

— v y> i T_ v y< i+1 to 
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y z,j+l/2 



V z ,i + V Z A+1 



B 



y,i+l/2 



B 



z,i+l/2 





2 




+ By,i+1 




2 


B z ,i 


+ B Z) i + i 
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(2-57) 
(2-58) 
(2-59) 



The state vector q at the cell center is updated by calculating the modified fluxes f x at the 
cell boundaries as follows: 



T n n+1 

1 



Qi 



At n 
~~Kx 



(fx,i-l/2 — /a;,i+l/ 2 )' 



x,i+l/2 



[FM)+FM+l)\ 



Ax 
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n ,1+1/2^,1+1/2' 



k=l 



( At n \ 

Pk,i+l/2 = Qk [ ~~^: a k,i+l/2 + 7fc,i+l/2 j a fc,i+l/2 ~ (ffc 
«fc,i+l/2 = ^fc,i+l/2 " (<7i+l ~~ 

r gM+i^ for ak , +i/2 + 0) 



9k,i+lJ 



7fc,i+l/2 



a fc,i+l/2 

for a fe)i+ i/ 2 = 0, 



# feji = sign(^ fc)i+ i/ 2 )max[0 3 inin{|^ )i4 .i/2|,5fc ) i-i/2sigii(5fc,t+i/2)}] 



ffjfc,i+l/2 



/At r 



V Ax 



a 



fc,i+l/2 



/At r 



I Ax «*,i+l/2 



a fc,i+l/2i 



(2-60) 

(2-61) 

(2-62) 
(2-63) 

(2-64) 
(2-65) 
(2-66) 



Qk(x) 



(2-67) 



^ + e fc for |x| < 2e k , 
\x\ for |x| > 2e fc . 

Since the use of contact steepener and rotational steepener produces spurious numerical oscillations 
in the adiabatic MHD code (RJ; RJF), we do not include the rotational steepener in the 
IMHD code. The time step size At n is restricted by the usual Courant condition for stability, 
At™ = C cour Ax/Max(|7;™ . +1/2 | + c£ i+1/2 ) with C cour < 1. 



2.3. Multi-dimensional Code 

We extend the one-dimensional IMHD code to more dimensions by using a Strang-type 
directional splitting (Strang 1968). Here we explain, as an example, the implementation of it in 
two-dimensional plane-parallel geometry. The two-dimensional IMHD equations written in the 
conservative form (Eq. |2-5| ) can be split into 

L x (x-sweep) : q t + ^ = 0, (2-68) 
L y (y-sweep) : q t + = 0. (2-69) 
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In a time step, we update the state vector q(x, y) along the x-direction with y fixed, followed along 
the y-direction with x fixed, 

q n+1 = L y L x q n . (2-70) 

In order to maintain a second-order accuracy in time, the order of directional sweeps is permuted 
in the next time step by L x L y . The time step size, At, is calculated at the start of the one 
complete sequence of L x L y L y L x and fixed through the sequence. 

In multi-dimensional simulations, numerical solutions may not satisfy V ■ B = due to 
discretization errors. Brackbill & Barnes (1980) pointed out that errors of non-zero V--B appear as 
a force parallel to the field. Non-zero V • B can be removed, for instances, either by incorporating 
an explicit divergence cleaning method as described in RJF or by implementing a scheme similar 
to the constrained transport scheme (Evans & Hawley 1988) which was described in details for the 
adiabatic MHD code in Ryu et al. (1998). Tests in the next section have been done using the 
explicit divergence cleaning method, and the next two paragraphs describe it briefly. 

At the beginning of MHD simulations, V • B = is satisfied. The updated magnetic field B, 
which is not in general divergence- free, can be decomposed as into two parts, 

B = -V<f> + V x V, (2-71) 

where <j) and V are scalar and vector functions respectively. Then the corrected magnetic field 
defined as B c = B + V</> becomes divergence-free. So the problem of the divergence-cleaning is 
reduced to find <p, which is described by the Poisson equation 

V 2 = —V • B. (2-72) 

In two-dimensional Cartesian geometry, for instance, the following finite difference representations 

B c _ d . , gj±lJ ~ gizlJ (2-7^ 



''U...J 'V'-./ • 2A^ . i 2 " 74 ) 



together with 



4>i+2,j — 2<t>i,j + <t>i-2,j _|_ 4>i,j+2 ~ + 4>i,j-2 _ _( B x ^ + ij — B x ^_ij ^ By^j+i — Byjj-i 



(2Ax) 2 (2Ay) 2 V 2Ax 2Ay 

(2-75) 

ensures 

DC DC DC DC 

2 Ax + 2 Ay " U ' l ^ bj 

within machine round-off error. Extensions to the three-dimension and/or to other geometry are 
straightforward. 

Eq. ( 2-75| ) is solved with boundary conditions specific to problems. In the problems with 
periodic boundaries (as in the decay test of the Alfven wave in §3.2), a fast Poisson solver based 
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on the fast Fourier transform can be used. In the nonlinear simulation of the Parker instability in 
§3.3 with reflection boundaries along one direction, the computational domain is doubled to that 
direction and the resulting boundaries are enforced to be periodic. In the two-dimensional shock 
tube tests in §3.1.2, doubling the computational domain in both directions also makes the resulting 
boundaries periodic. Note that in Eq. (2-75) 0's are coupled with those at every other cell in a 
column and row. So the (extended) computational domain is divided into four sub-domains, and 
0's are computed in those sub-domains separately. 



3. TESTS 

In this section we present the results of three tests. The first and the second are isothermal 
versions of MHD shock tubes and decay of an Alfven wave, respectively (RJ; RJF). The shock 
tube test shows the ability of the IMHD code to handle all the three MHD wave family structures, 
while the decay test of an Alfven wave measures numerical dissipation. The third test is the 
simulation of a real astrophysical situation, the nonlinear evolution of the Parker instability under 
a uniform gravity. In all the tests, we set the isothermal speed a = 1. 



3.1. Shock Tube Tests 

Based on the work of RJ, we have devised four shock tube problems which include 
discontinuities and rarefaction waves of IMHDs. To confirm the validity of our numerical solutions 
we have compared them to the analytic solutions obtained with an exact, nonlinear MHD Riemann 
solver described in RJ. That Riemann solver iterates from an initial guess of the solution for 
the full set of MHD waves based on the given left and right states. Iteration continues until the 
solutions to the innermost wave zone reached from the two opposite directions agree within some 
specified limit (in practice, a relative error 10~ 5 ). The four shock tube solutions that we have 
applied in the test described below are listed in Table 1. C COU r = 0.8 and ei = 0.3 (for fast modes), 
€2 = 0.0 (for Alfven modes), e% = 0.3 (for slow modes) have been used in the shock tube test 
calculations. 



3.1.1. One- Dimensional Shock Tube Tests 

The one-dimensional simulations of the shock tube problems have been done with 512 cells in 
a computational tube bounded by x = [0, 1]. We plot in following figures the resulting p, B y , B z , 
v x , v y , and v z at each cell with open circles and the analytic solutions with lines. 

Figure la shows the result of the first shock tube problem at t = 0.1 with the initial 
condition of a left state (p = l,v x = 0,v y = 0, v z = 0,B y = 5/\/4vr, B z = 0), a right state 
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(p = 0.1, v x = 0, v y = 0,v z = 0, By = 2/\/47r, B z = 0), and B x = 3/\/47r. It exhibits the 
capturing of a fast rarefaction wave, a slow rarefaction wave, a slow shock, and a fast shock 
whose structures are plotted in the figure from left to right. There is no contact discontinuity. 
The fast and slow shocks are resolved sharply within several cells. In order to see the capturing 
rotational discontinuities, we have set up the initial condition of the second shock tube problem 
as: a left state (p = 1.08, v x = 1.2, v y = 0.01, v z = 0.5, B y = 3.6/\/47r, B z = 2/v / 47r), a right 
state (p = l,v x = 0,v y = 0,v z = 0,B y = 4/y/4n,B z = 2/y/4n), and B x = 2/y/4n. Figure lb 
shows the result at t = 0.2. There are two fast shocks propagating outmost, and two slow 
shocks interior to those. Two rotational discontinuities lie between the fast and slow shocks. 
Here the strong fast shocks are resolved sharply, but the weak slow shocks and rotational 
discontinuities spread over more cells. The third shock tube problem has been set up with 
the initial condition of a left state (p = 0.12, v x = 24, v y = 0,v z = 0,B y = 3/V4tt,B z = 0), 
a right state (p = 0.3, v x = —15,v y = 0, v z = 0,B y = 0,B Z = 3/\/47r), and B x = 0. Two 
oppositely moving magnetosonic shocks and a tangential discontinuity at t = 0.2 are shown in 
Figure lc. The magnetosonic shocks are again resolved sharply, but the tangential discontinuity 
spreads over ~ 20 cells. In the fourth shock tube problem, the initial condition has been 
set up with a left state (p = l,v x = —l,v y = 0,v z = 0,B y = 1,B Z = 0), a right state 
(p = l,v x = l,v y = 0, v z = 0,B y = 1,B Z = 0), and B x = 0. The result at t = 0.16 in Figure Id 
shows two oppositely-moving identical magnetosonic rarefactions. 

3.1.2. Two-Dimensional Shock Tube Tests 

The two-dimensional simulations of the shock tube problems have been done with 256 x 256 
cells in a computational domain bounded by x = [0, 1] and y = [0, 1]. Initially, the domain is 
divided into two parts by a diagonal line joining the two points (0,1) and (1,0). The left state 
of the initial conditions for the one-dimensional shock tube problems has been assigned to the 
lower left part and the right state to the upper right part. The generated structures, including 
discontinuities and rarefactions, propagate parallel to the other diagonal line joining the two 
points (0,0) and (1,1). 

In Figures 2a and 2b, two-dimensional correspondences of Figures la and lb are plotted. 
In the figures, the following subscripts are used: || for parallel components of velocity and 
magnetic field along the diagonal line joining the two points (0,0) and (1,1), -L for perpendicular 
components but still in the computational plane, and z for components out of the plane. Although 
the resolution of the two-dimensional simulations, 256 x 256 cells, is lower than that of the 
one-dimensional ones, 512 cells, we see all the structures have been captured correctly. 
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3.2. Decay of an Alfven Wave 



RJF carried out a test of the decay of linear waves in order to estimate numerical dissipations 
(resistivity and viscosity) in their adiabatic MHD code. Following the same idea, the decay of 
a linear Alfven wave has been calculated and numerical dissipation in our IMHD code has been 
estimated. The IMHD equations for viscous and resistive fluid can be written as 



! + v. W = o, 

dv 1 1 1 

— + v ■ Vv + -V(a 2 p) - -(V x B) x B = -d k a ik , 

at p p p 



dB 
~dt 



V x (v x B) = r]V 2 B. 



In the momentum equation, the viscosity tensor is given by, 

Oik = p{dkVi + div k - p ik V • v) + Q5 ik V ■ v, 



(3-1) 
(3-2) 
(3-3) 

(3-4) 



where p and ( are the dynamic shear and bulk viscosity, and r\ is the electrical resistivity. Under 
uniform density, po, and uniform magnetic field, B = Bqx, the complex angular frequency of 
Alfven waves is predicted from the linear analysis to be 



j = - [ — + ri) k z ± cj^k 
Po 



1 



1 



p 

— v 

Po 



.1/2 



(3-5) 



where ca = \j Bq/2pq is the Alfven speed along the wave propagation direction and k = (k^ + ky) 1 ^ 2 
is the total wavenumber. Note that the complex angular frequency of the isothermal Alfven waves 
is exactly the same as that of the adiabatic ones (RJF). We define the decay rate as 



r A = 



l 



p_ 

Po 



+ t) \k 



(3-6) 



For the decay test of a linear Alfven wave with the IMHD code, we have set up an initial 
condition such that, po = 1> $ v z = ^amp c A sin(k x x + k y y), B = 1 • x, and all other quantities are 
equal to zero. The calculations have been done in a square periodic box with size L = 1 using 
from 8x8 cells to 128 x 128 cells by increasing twice the number of cells in each direction. We set 
k x = k y = 2it/L. Numerical parameters used are ei = £2 = £3 = and C COU r = 0.9 (the result is not 
sensitive to these values). Figure 3a shows the decay of the Alfven wave calculated with 32 x 32 
cells. By fitting the peak points of the decay pattern with respect to time, we have estimated 
decay rate. In Figure 3b the resulting normalized decay rates as well as Reynolds numbers (see 
RJF for definition) are shown. Numerical Reynolds numbers scale almost as R oc n^ ell indicating 
the code has a second-order accuracy. Compared to the adiabatic MHD code, the IMHD code has 
smaller (up to 50%) numerical dissipation. This is partly because the IMHD code has one less 
mode (entropy mode). 
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3.3. Parker Instability under a Uniform Gravity 

Nonlinear development of the Parker instability under a point-mass dominated gravity was 
simulated by Matsumoto and his collaborators (Matsumoto et al. 1988; Matsumoto et al. 1990; 
Matsumoto & Shibata 1992). And recently, Basu et al. (1996, 1997) simulated the nonlinear 
evolution of the Parker instability under a uniform gravity. As the final test of our IMHD code, 
we have also followed the nonlinear evolution of the Parker instability under the uniform gravity. 
By comparing our results with those in Basu et al. (1996), the code's ability to handle a practical 
problem of astrophysics can be proved. 

Since the Parker system is assumed initially to be in an isothermal equilibrium, an IMHD 
code is a natural tool for simulations. In the IMHD equations the externally given gravity -gz is 
treated as a source term and placed on the right hand side of Eq. (|2-5[ ) with the source vector 
defined by S = (0,0,0, — gv z , 0, 0, 0) T . The gravity has the z-component only, so the source term 
is evaluated only when the state vector is updated along the z-axis 

Since we use the Strang-type directional splitting in order to reduce multi-dimensional problems 
to one-dimensional ones, we also use the same technique to split the hyperbolic system with a 
source term into two parts, 

Part A: ^ + ^ = 0, (3-8) 
at oz 

Part B: ^ = S. (3-9) 
dt v ; 

Part A is solved by the TVD algorithm, and Part B by a forward-time-difference. To minimize a 
'splitting error', Part A and Part B are solved by a BAB sequence with time step size 0.5At for 
Part B and At for Part A. The step size At is determined from the Courant condition. 

The Parker system composed of isothermal gas and magnetic field yB^z) takes under the 
uniform gravity an equilibrium configuration given by 

where the gas scale height is H = (1 + a)a 2 /g and the initial ratio of the magnetic to the gas 
pressure a(= i?g/[2poo 2 ]) is assumed a constant. We have chosen a = 1.0 in the simulation. 

The computation domain covers < y < 12H and < z < 12H. According to the linear 
stability analysis \2H is the horizontal wavelength corresponding to the maximum growth rate 
(Parker 1966). Periodic condition has been used in the y-boundaries, while reflecting condition in 
the z-boundaries. The density scale height H , the isothermal sound speed a, the initial midplane 
density po(0), and the initial midplane field strength B o (0) have been chosen as the units of length, 
velocity, density, and magnetic field, respectively. 
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To initiate the instability we have added random velocity perturbations to the equilibrium 
profile of Eq. ( |3-10 ). Standard deviation of the perturbation velocity is 10 _4 a for each of the 



velocity components. To check whether the system follows, in the initial stage, the prescription of 
the linear analysis, we plot the logarithmic values of the root-mean-square velocity against time. 
In Figure 4a the dotted and dashed lines are for the horizontal and vertical components of the 
velocity, respectively. In the same figure the solid line represents the linear growth of a rate 0.34, 
which is the maximum growth rate of the system. At the very early stage the system undergoes 
a transient phase of adjustment, and then quickly develops the Parker instability at the rate 
predicted by the linear analysis. The linear growth gets saturated near t ~ 40. 

The whole development of the Parker instability may be divided into three phases: The linear 
phase lasts up to t ~40, from then on the system undergoes the nonlinear phase until t ~ 57, 
and finally it reaches the damping oscillatory phase. Iso-contours and grey maps for density (left 
panels) and magnetic field lines and velocity vectors (right panels), in Figure 4b, present the snap 
shots of the system at the end of the linear phase (t = 40), at the end of the nonlinear phase 
(t = 57), and finally at t = 80 of the damping oscillatory phase. 

In the linear phase the perturbations grow predominantly in the upper region. In the nonlinear 
phase the perturbations gradually move towards the midplane. Through the linear and nonlinear 
phases, our simulation renders features that closely agree with those of Basu et al. (1996). As more 
matter accumulates, already compressed gas in the valley gets over-compressed. The increased gas 
pressure bounces the valley matter back to the upper region, and at the same time, the built-up 
pressure at the valley gets somewhat eased off. This in turn brings the matter back to the valley. 
The system now enters the oscillatory phase of the Parker instability. As the field lines are pushed 
deep down to the valley by the weight of over-lying matter, the curvature of the lines becomes 
small to the degree that magnetic field lines undergo reconnection. Due to the reconnection the 
matter drops down off the reconnected line, thereby the matter is allowed to move across the 
magnetic field line. The field line is now relieved from the burden of supporting the gas against 
the external gravity, and floats upwards. On the other hand, the field line located just below 
the reconnected one has to support more weight than before. Consequently this line now gets 
reconnected. In this way a redistribution of matter with respect to the field lines continues to 
occur, until there is no more reconnection. The system finally settles in an equilibrium. Since the 
reconnection drives the system to violate the flux-freezing condition, the final configuration of the 
system is different from that of the Mouschovias equilibrium (1974). 



4. CONCLUSIONS 

We have developed one- and multi-dimensional numerical codes to solve the IMHD equations, 
which are isothermal analogues of the previous adiabatic codes (RJ; RJF). Both the isothermal 
and adiabatic codes are based on the same scheme, an explicit finite-difference scheme on an 
Eulerian grid called TVD, which is a second-order-accurate extension of the Roe- type upwind 
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scheme. The shock tube tests have showed that both codes capture correctly all the structures in 
MHDs. From the decay test of a linear Alfven wave, we have found that numerical dissipation of 
the IMHD code is somewhat smaller than that of the adiabatic MHD code. 

The robustness of the adiabatic code has been demonstrated through the simulations of 
MHD flows such as the Kelvin-Helmholtz instability (Frank et al. 1996; Jones et al. 1997) and 
jets (Frank et al. 1998; Jones et al. 1998), and that of the isothermal code has been done through 
the simulation of the Parker instability under the uniform gravity in this paper. Furthermore, 
both codes are fast enough to simulate multi-dimensional, astrophysical MHD flows using modest 
computational resources. Both codes run at about 400 MFlops on a Cray C90 processor (RJF), 
and the isothermal code updates zones about twice as fast as the adiabatic code. Together 
with the adiabatic code, the isothermal code is a useful tool to study the nonlinear evolution of 
astrophysical MHD flows. 
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FIGURES 

Fig. la. — One-dimensional IMHD shock tube test. The initial condition is (p, v x ,v y , v z ,B y , B z ) = 
(l,0,0,0,5/\/4vf,0) in the left region, (p, v x , v y , v z , B y , B z ) = (0.1, 0, 0, 0, 2/a/4tt, 0) in the right 
region, B x = 2>/\f4m and a = 1 for the whole computational interval. Open circles represent the 
numerical solution, while lines represent the analytic solution with an exact nonlinear Riemann 
solver. The calculation has been done with 512 cells. A snapshot at t = 0.1 shows from left to right 
(1) fast rarefaction, (2) slow rarefaction, (3) slow shock, and (4) fast shock. 

Fig. lb. — One-dimensional IMHD shock tube 

test. The initial condition is (p,v x ,v y ,v z , B y , B z ) = (1.08, 1.2, 0.01, 0.5, 3.6/\/47r, 2/\/47r) in the 
left region, (p, v x , v y , v z , B y , B z ) = (1,0,0,0,4/^/4^,2/^/47?) in the right region, B x = 2/^/4tt and 
a = 1 for the whole computational interval. Open circles represent the numerical solution, while 
lines represent the analytic solution with an exact nonlinear Riemann solver. The calculation has 
been done with 512 cells. A snapshot at t = 0.2 shows from left to right (1) fast shock, (2) rotational 
discontinuity, (3) slow shock, (4) slow shock, (5) rotational discontinuity, and (6) fast shock. 

Fig. lc. — One-dimensional IMHD shock tube test. The initial condition is (p, v x ,v y , v z ,B y , B z ) = 
(0.12,24,0,0,3/^/4^,0) in the left region, (p, v x , v y , v z , B y , B z ) = (0.3,-15,0,0,0,3/^/4^) in the 
right region, B x = and a = 1 for the whole computational interval. Open circles represent the 
numerical solution, while lines represent the analytic solution with an exact nonlinear Riemann 
solver. The calculation has been done with 512 cells. A snapshot at t = 0.2 shows from left to right 
(1) magnetosonic shock, (2) tangential discontinuity, and (3) magnetosonic shock. 

Fig. Id. — One-dimensional IMHD shock tube test. The initial condition is (p, v x ,v y , v z ,B y , B z ) = 
(1, —1,0, 0, 1, 0) in the left region, (p, v x ,v y ,v z , B y , B z ) = (1, 1, 0, 0, 1, 0) in the right region, B x = 
and a = 1 for the whole computational interval. Open circles represent the numerical solution, 
while lines represent the analytic solution with an exact nonlinear Riemann solver. The calculation 
has been done with 512 cells. A snapshot at t = 0.16 shows from left to right (1) magnetosonic 
rarefaction, and (2) magnetosonic rarefaction. 

Fig. 2a. — Two-dimensional IMHD shock tube test. The initial condition is (p, vy, v±, v z , B±, B z ) = 
(1,0,0,0,5/V4tt,0) in the lower left region, (p, v\\, v±, v z , B ± , B z ) = (0.1, 0, 0, 0, 2/ v / 4vf, 0) in the 
upper right region, Bii = 3/\/47r and a = 1 in the whole computational domain. Open circles 
represent the numerical solution, while lines represent the analytic solution with an exact nonlinear 
Riemann solver. The calculation has been done with 256 x 256 cells. The structures shown at 
t = 0.1\/2 along a diagonal line joining the two points (0,0) to (1,1) are same as those of Figure la. 

Fig. 2b. — Two-dimensional IMHD shock tube test. The 

initial condition is (p, v\\ , v±, v z , B±, B z ) = (1.08, 1.2, 0.01, 0.5, 3.6/\/47r, 2/\/47r) in the lower left 
region, (p,v\\,v±,v z ,B±,B z ) = (1, 0, 0, 0, A/y/Air, 2/\f4m) in the upper right region, B\\ = 2/\f4m 
and a = 1 in the whole computational domain. Open circles represent the numerical solution, while 
lines represent the analytic solution with an exact nonlinear Riemann solver. The calculation has 
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been done with 256 x 256 cells. The structures shown at t = 0.2a/2 along a diagonal line joining 
the two points (0,0) to (1,1) are same as those of Figure lb. 

Fig. 3a. — Time evolution of < SB 2 > x l 2 and < 5v 2 > x l 2 in the decay test of a two-dimensional 
Alfven wave. Initially, a standing Alfven wave has been set up in a computational domain with 
32 x 32 cells, and its oscillation has been followed. 

Fig. 3b. — Normalized decay rate, TaL/ca, and magnetic Reynolds number, R, as a function 
of the number of cells along one direction of the computation domain. At a given resolution, the 
peak-to-peak decay rate of the root-mean-square of z-velocity (top) and the corresponding Reynolds 
number (bottom) are plotted with filled circles, respectively. The calculations have been done with 
8x8, 16 x 16, 32 x 32, 64 x 64, and 128 x 128 cells. For comparison, dotted lines of (TaL/ca) oc n~ 2 xl 
and R oc n 2 cU are drawn. 

Fig. 4a. — Time evolution of the root-mean-square of the horizontal velocity, < v 2 > x / 2 , and the 
vertical velocity, < v 2 > 1 / 2 , in a simulation of the Parker instability under a uniform gravity. The 
magnetohydrostatic equilibrium state together with random velocity perturbations has been given 
as an initial condition of the simulation in the computational domain of 256 x 256 cells. The solid 
line represents the predicted maximum linear growth with perturbation wavelength \ y = 12, and 
X z /2 = 12. The normalization units are the isothermal sound speed and the scale height. 

Fig. 4b. — Evolution of the Parker instability under a uniform gravity. At three time epochs t = 40 
(top), t = 57 (middle), and t = 80 (bottom), grey maps of density together with equi-density lines 
are plotted in left panels, and the velocity vectors with magnetic field lines in right panels. The 
values of the ten equi-density lines are the initial exponential densities at z = 1, ■ ■ ■ , 10. Magnetic 
field lines are chosen so that the magnetic flux between two consecutive lines is constant. At each 
time epoch, the unit of the velocity vectors is represented. 
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Table la. Shock Tube Test la 



p 


Vx 


Vy 




By 


B z 


1.0000E+00 


0.0000E+00 


0.0000E+00 


0.0000E+00 


1.4105E+00 


0.0000E+00 


5.7648E-01 


9.3200E-01 


-5.3737E-01 


0.0000E+00 


5.9825E-01 


0.0000E+00 


3.0968E-01 


1.3718E+00 


-1.0767E-02 


0.0000E+00 


7.8902E-01 


0.0000E+00 


1.2358E-01 


7.2565E-01 


-7.6338E-01 


0.0000E+00 


9.0720E-01 


0.0000E+00 


1.0000E-01 


0.0000E+00 


0.0000E+00 


0.0000E+00 


5.6419E-01 


0.0000E+00 



with B x = 8.4628S - 01 



Table lb. Shock Tube Test lb 



p 


Vx 


Vy 


Vz 


By 


B z 


1.0800E+00 


1.2000E+00 


1.0000E-02 


5.0000E-01 


1.0155E+00 


5.6419E-01 


1.5087E+00 


6.4673E-01 


1.3132E-01 


5.6740E-01 


1.4677E+00 


8.1542E-01 


1.5087E+00 


6.4673E-01 


2.4196E-01 


3.0857E-01 


1.6036E+00 


4.9750E-01 


1.7451E+00 


6.0765E-01 


7.3388E-02 


2.5628E-01 


1.4736E+00 


4.5716E-01 


1.3560E+00 


5.4030E-01 


-2.1440E-01 


1.6699E-01 


1.6825E+00 


5.2198E-01 


1.3560E+00 


5.4030E-01 


-1.2262E-01 


-6.1311E-02 


1.5757E+00 


7.8783E-01 


1.0000E+00 


0.0000E+00 


0.0000E+00 


0.0000E+00 


1.1284E+00 


5.6419E-01 



with B x = 5.6419S - 01 
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Table lc. Shock Tube Test lc 



p 


Vx 


Vy 


Vz 


By 


B z 


1.2000E-01 
1.7079E+00 
4.1960E+00 
3.0000E-01 


2.4000E+01 
9.2149E-02 
9.2149E-02 

-1.5000E+01 


2.3130E-16 
2.3130E-16 
0.0000E+00 
0.0000E+00 


0.0000E+00 
0.0000E+00 
0.0000E+00 
0.0000E+00 


8.4628E-01 
1.2045E+01 
7.2478E-16 
0.0000E+00 


4.4409E-16 
6.3206E-15 
1.1837E+01 
8.4628E-01 



with B x = 0.0000S + 00 



Table Id. Shock Tube Test Id 



P V X Vy V Z By B Z 



1.0000E+00 -1.0000E+00 0.0000E+00 0.0000E+00 1.0000E+00 0.0000E+00 

4.6392E-01 7.8159E-08 0.0000E+00 0.0000E+00 4.6392E-01 0.0000E+00 

1.0000E+00 1.0000E+00 0.0000E+00 0.0000E+00 1.0000E+00 0.0000E+00 



with B x = 0.0000S + 00 
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